Skewness Test (D’Agostino)

The skewness test checks whether a sample shows statistically significant asymmetry.

It is commonly used as a diagnostic for:

• symmetry assumptions (many parametric procedures implicitly assume roughly symmetric errors)

• whether a transformation (log/sqrt/Box–Cox) might be useful

• as one component of normality diagnostics (but it is not a full normality test by itself)

Learning goals

By the end you should be able to:

• explain what skewness measures (right/left tail intuition)

• state the hypotheses behind the skewness test

• implement D’Agostino’s skewness test with NumPy only

• interpret the test statistic and p-value (including one-sided alternatives)

• understand common pitfalls (outliers, sample size, “significant but tiny” effects)

import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=6, suppress=True)

rng = np.random.default_rng(42)

Prerequisites

• Sample mean and central moments

• Z-scores and p-values

• The standard normal distribution

If any of these are fuzzy, you can still follow the notebook: the code is written to make each step explicit.

1) Intuition: what skewness measures

Skewness is a signed measure of asymmetry:

• Positive skewness → a longer right tail (“right-skewed”)

• Negative skewness → a longer left tail (“left-skewed”)

A common (moment-based) sample skewness is

[ ext{skew}(x) = rac{m_3}{m_2^{3/2}}, \qquad m_k = rac{1}{n}\sum_{i=1}^n (x_i - ar{x})^k. ]

This uses the third central moment (m_3), so it is very sensitive to outliers.

def sample_skewness(x: np.ndarray) -> float:
    "Moment-based sample skewness: m3 / m2^(3/2) using 1/n moments."
    x = np.asarray(x, dtype=float).ravel()
    x = x[np.isfinite(x)]
    n = x.size
    if n < 3:
        return float("nan")

    mean = x.mean()
    centered = x - mean
    m2 = np.mean(centered**2)
    if m2 == 0:
        return 0.0
    m3 = np.mean(centered**3)
    return float(m3 / (m2 ** 1.5))


n_demo = 4000
x_sym = rng.normal(0, 1, size=n_demo)
x_right = rng.exponential(scale=1.0, size=n_demo)
x_left = -rng.exponential(scale=1.0, size=n_demo)

skews = {
    "Symmetric (normal)": sample_skewness(x_sym),
    "Right-skewed (exponential)": sample_skewness(x_right),
    "Left-skewed (-exponential)": sample_skewness(x_left),
}

skews

{'Symmetric (normal)': 0.003983194961544123,
 'Right-skewed (exponential)': 1.979806008441165,
 'Left-skewed (-exponential)': -1.901827954731168}

fig = make_subplots(
    rows=1,
    cols=3,
    subplot_titles=[f"{k}<br>skew={v:.3f}" for k, v in skews.items()],
)

datasets = [
    ("Symmetric (normal)", x_sym),
    ("Right-skewed (exponential)", x_right),
    ("Left-skewed (-exponential)", x_left),
]

for col, (name, x) in enumerate(datasets, start=1):
    fig.add_trace(
        go.Histogram(
            x=x,
            nbinsx=80,
            histnorm="probability density",
            name=name,
            marker=dict(line=dict(width=0)),
        ),
        row=1,
        col=col,
    )

fig.update_layout(
    height=320,
    width=1000,
    showlegend=False,
    title_text="Skewness: symmetric vs right/left skew",
)
fig.update_xaxes(range=[-5, 5])
fig.show()

2) What the skewness test is (and what it is not)

Goal

We want to know whether the skewness we observe could plausibly be due to sampling noise when the underlying distribution is normal (which has population skewness 0).

Hypotheses

For the standard two-sided skewness test:

• (H_0): the population skewness is the same as the normal distribution (i.e. 0)

• (H_1): the population skewness is not 0

You can also run one-sided versions:

• alternative='greater': (H_1): skewness (> 0) (right-skew)

• alternative='less': (H_1): skewness (< 0) (left-skew)

What it is not

• It is not a complete normality test. A distribution can be symmetric (skew ≈ 0) but still very non-normal (e.g. heavy-tailed).

• It does not tell you why data are skewed (mixtures, censoring, outliers, bounded support, etc.). You should always visualize the data.

3) D’Agostino’s skewness test: the core idea

D’Agostino’s skewness test takes the sample skewness (b_2) and transforms it into a value (Z) that is approximately standard normal under (H_0):

[ Z pprox \mathcal{N}(0,1) \quad ext{(when the data are normal and } n\ge 8 ext{)}. ]

That lets us compute a p-value just like a z-test.

The transformation (as used by SciPy’s skewtest)

Let (n) be the sample size and (b_2) the (moment-based) sample skewness. Define

[ y = b_2\sqrt{rac{(n+1)(n+3)}{6(n-2)}}. ]

Then compute

[ eta_2 = rac{3,(n^2 + 27n - 70)(n+1)(n+3)}{(n-2)(n+5)(n+7)(n+9)}, \quad W^2 = -1 + \sqrt{2(eta_2 - 1)}. ]

and

[ \delta = rac{1}{\sqrt{ frac{1}{2}\ln(W^2)}}, \quad lpha = \sqrt{rac{2}{W^2 - 1}}. ]

Finally,

[ Z = \delta,\operatorname{asinh}!\left(rac{y}{lpha} ight) = \delta,\ln\left(rac{y}{lpha} + \sqrt{\left(rac{y}{lpha} ight)^2 + 1} ight). ]

Why the sample size constraint?

The approximation is derived under normality and is intended for (n\ge 8). For smaller samples, you should prefer resampling-based checks or simply treat skewness as an effect size.

4) NumPy-only implementation

The functions below implement:

• sample skewness (b_2)

• D’Agostino’s transformed statistic (Z)

• p-values using a NumPy approximation to the normal CDF (via an erf approximation)

The goal here is transparency, not squeezing out every last bit of numerical accuracy.

def erf_approx(x: np.ndarray) -> np.ndarray:
    "Vectorized erf approximation (Abramowitz & Stegun 7.1.26)."
    x = np.asarray(x, dtype=float)

    sign = np.sign(x)
    ax = np.abs(x)

    p = 0.3275911
    t = 1.0 / (1.0 + p * ax)

    a1 = 0.254829592
    a2 = -0.284496736
    a3 = 1.421413741
    a4 = -1.453152027
    a5 = 1.061405429

    poly = (((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t)
    y = 1.0 - poly * np.exp(-ax * ax)
    return sign * y


def norm_cdf(z: np.ndarray) -> np.ndarray:
    "Standard normal CDF Φ(z) using erf approximation."
    z = np.asarray(z, dtype=float)
    return 0.5 * (1.0 + erf_approx(z / np.sqrt(2.0)))


def pvalue_from_z(z: np.ndarray, alternative: str = "two-sided") -> np.ndarray:
    "Convert a z-statistic to a p-value under N(0,1)."
    z = np.asarray(z, dtype=float)

    if alternative == "two-sided":
        p = 2.0 * norm_cdf(-np.abs(z))
    elif alternative == "greater":
        p = norm_cdf(-z)
    elif alternative == "less":
        p = norm_cdf(z)
    else:
        raise ValueError("alternative must be one of: 'two-sided', 'greater', 'less'")

    return np.clip(p, 0.0, 1.0)


def dagostino_z_from_skewness(b2: np.ndarray, n: int) -> np.ndarray:
    "D'Agostino transform from sample skewness b2 to Z (approx N(0,1) under H0)."
    if n < 8:
        raise ValueError("D'Agostino skewness test requires n >= 8")

    b2 = np.asarray(b2, dtype=float)
    nf = float(n)

    y = b2 * np.sqrt(((nf + 1.0) * (nf + 3.0)) / (6.0 * (nf - 2.0)))

    beta2 = (
        3.0
        * (nf**2 + 27.0 * nf - 70.0)
        * (nf + 1.0)
        * (nf + 3.0)
        / ((nf - 2.0) * (nf + 5.0) * (nf + 7.0) * (nf + 9.0))
    )

    w2 = -1.0 + np.sqrt(2.0 * (beta2 - 1.0))
    delta = 1.0 / np.sqrt(0.5 * np.log(w2))
    alpha = np.sqrt(2.0 / (w2 - 1.0))

    # asinh(u) == log(u + sqrt(u^2 + 1)) but is more numerically stable
    return delta * np.arcsinh(y / alpha)


def dagostino_skewtest(x: np.ndarray, alternative: str = "two-sided") -> dict:
    "D'Agostino skewness test (NumPy-only implementation)."
    x = np.asarray(x, dtype=float).ravel()
    x = x[np.isfinite(x)]
    n = x.size
    if n < 8:
        raise ValueError("D'Agostino skewness test requires at least 8 observations")

    b2 = sample_skewness(x)
    z = dagostino_z_from_skewness(b2, n)
    p = pvalue_from_z(z, alternative=alternative)

    return {
        "n": int(n),
        "skewness": float(b2),
        "z": float(np.asarray(z)),
        "pvalue": float(np.asarray(p)),
        "alternative": alternative,
    }

5) Worked examples + interpretation

Interpretation guide (two-sided):

• Small p-value (e.g. < 0.05) → evidence the data are asymmetric (skewness ≠ 0) relative to normal sampling noise.

• The sign of Z tells you the direction:

  • Z > 0 → right-skew (long right tail)

  • Z < 0 → left-skew (long left tail)

Remember: with large (n), even tiny skewness can become “statistically significant”. Always look at the effect size (the skewness value) and the plot.

def summarize_test(x: np.ndarray, name: str, alternative: str = "two-sided") -> dict:
    res = dagostino_skewtest(x, alternative=alternative)
    return {
        "dataset": name,
        "n": res["n"],
        "skewness": res["skewness"],
        "z": res["z"],
        "pvalue": res["pvalue"],
        "alternative": res["alternative"],
    }


n = 60
x1 = rng.normal(0, 1, size=n)
x2 = rng.exponential(1.0, size=n)
x3 = -rng.exponential(1.0, size=n)

results = [
    summarize_test(x1, "Normal sample"),
    summarize_test(x2, "Right-skew (exponential)"),
    summarize_test(x3, "Left-skew (-exponential)"),
]

results

[{'dataset': 'Normal sample',
  'n': 60,
  'skewness': 0.33403897478489747,
  'z': 1.1391470464804974,
  'pvalue': 0.25464193619392184,
  'alternative': 'two-sided'},
 {'dataset': 'Right-skew (exponential)',
  'n': 60,
  'skewness': 1.0681651963774472,
  'z': 3.204306467065742,
  'pvalue': 0.001354010697767749,
  'alternative': 'two-sided'},
 {'dataset': 'Left-skew (-exponential)',
  'n': 60,
  'skewness': -2.314018629424382,
  'z': -5.399094653879995,
  'pvalue': 6.713003297686981e-08,
  'alternative': 'two-sided'}]

fig = make_subplots(
    rows=1,
    cols=3,
    subplot_titles=[
        f"{r['dataset']}<br>skew={r['skewness']:.3f}, Z={r['z']:.2f}, p={r['pvalue']:.3g}" for r in results
    ],
)

for col, (x, r) in enumerate([(x1, results[0]), (x2, results[1]), (x3, results[2])], start=1):
    fig.add_trace(
        go.Histogram(
            x=x,
            nbinsx=30,
            histnorm="probability density",
            name=r["dataset"],
            marker=dict(line=dict(width=0)),
        ),
        row=1,
        col=col,
    )

fig.update_layout(
    height=320,
    width=1100,
    showlegend=False,
    title_text="Skewness test on small samples (n=60)",
)
fig.show()

One-sided alternatives

If you have a directional question, you can use a one-sided test:

• greater: “is it right-skewed?” (skewness > 0)

• less: “is it left-skewed?” (skewness < 0)

One-sided p-values are smaller only when the observed skew is in the hypothesized direction.

x = rng.exponential(1.0, size=60)

res_two = dagostino_skewtest(x, alternative="two-sided")
res_greater = dagostino_skewtest(x, alternative="greater")
res_less = dagostino_skewtest(x, alternative="less")

{
    "two-sided": res_two,
    "greater": res_greater,
    "less": res_less,
}

{'two-sided': {'n': 60,
  'skewness': 2.616686399904669,
  'z': 5.78336467838644,
  'pvalue': 7.344590824409636e-09,
  'alternative': 'two-sided'},
 'greater': {'n': 60,
  'skewness': 2.616686399904669,
  'z': 5.78336467838644,
  'pvalue': 3.672295412204818e-09,
  'alternative': 'greater'},
 'less': {'n': 60,
  'skewness': 2.616686399904669,
  'z': 5.78336467838644,
  'pvalue': 0.9999999963277046,
  'alternative': 'less'}}

6) Why sample size matters (p-value vs skewness)

Because the test statistic is a standardized measure of skewness, larger samples make it easier to detect small asymmetries.

The plot below shows the two-sided p-value as a function of sample skewness (b_2) for different (n).

b2_grid = np.linspace(-2.5, 2.5, 501)

fig = go.Figure()

for n in [10, 30, 100, 500]:
    z_grid = dagostino_z_from_skewness(b2_grid, n=n)
    p_grid = pvalue_from_z(z_grid, alternative="two-sided")
    fig.add_trace(go.Scatter(x=b2_grid, y=p_grid, mode="lines", name=f"n={n}"))

fig.add_hline(y=0.05, line_dash="dash", line_color="black")

fig.update_layout(
    title="Two-sided p-value vs sample skewness (D'Agostino transform)",
    xaxis_title="sample skewness b2",
    yaxis_title="p-value",
    yaxis_type="log",
    height=420,
    width=900,
)
fig.show()

7) Does Z really look standard normal under H0?

Under the null (normal data), the transformed statistic (Z) is designed to be approximately (\mathcal{N}(0,1)).

We can sanity-check that with a small Monte Carlo simulation.

def skewness_vectorized(samples: np.ndarray) -> np.ndarray:
    "Vectorized skewness for samples shaped (n_sims, n)."
    samples = np.asarray(samples, dtype=float)
    mean = samples.mean(axis=1, keepdims=True)
    centered = samples - mean
    m2 = np.mean(centered**2, axis=1)
    m3 = np.mean(centered**3, axis=1)

    out = np.zeros_like(m2)
    mask = m2 > 0
    out[mask] = m3[mask] / (m2[mask] ** 1.5)
    out[~mask] = 0.0
    return out


n = 20
n_sims = 10_000
samples = rng.normal(0, 1, size=(n_sims, n))

b2 = skewness_vectorized(samples)
z = dagostino_z_from_skewness(b2, n=n)

z_mean = float(np.mean(z))
z_std = float(np.std(z, ddof=0))
rejection_rate = float(np.mean(pvalue_from_z(z, alternative="two-sided") < 0.05))

{"z_mean": z_mean, "z_std": z_std, "rejection_rate@0.05": rejection_rate}

{'z_mean': 0.010093399552716806,
 'z_std': 0.9999562884320329,
 'rejection_rate@0.05': 0.0476}

# Histogram of Z with a standard normal PDF overlay

x_grid = np.linspace(-4, 4, 400)
pdf = (1.0 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * x_grid**2)

fig = go.Figure()
fig.add_trace(
    go.Histogram(
        x=z,
        nbinsx=60,
        histnorm="probability density",
        name="Simulated Z (H0: normal)",
        opacity=0.75,
        marker=dict(line=dict(width=0)),
    )
)
fig.add_trace(go.Scatter(x=x_grid, y=pdf, mode="lines", name="N(0,1) PDF", line=dict(color="black")))

fig.update_layout(
    title=f"Under H0, Z is approximately N(0,1) (n={n}, sims={n_sims})",
    xaxis_title="Z",
    yaxis_title="density",
    height=420,
    width=900,
)
fig.show()

8) Practical notes and pitfalls

• Outliers can dominate: skewness uses ((x-ar{x})^3). One extreme point can flip the conclusion.

• Statistical vs practical significance: with large (n), tiny skewness can yield tiny p-values.

• Not a full normality test: skewness ≈ 0 does not imply normality.

• Independence matters: if observations are dependent (time series, clustered data), p-values can be misleading.

• Use visuals: combine the test with histograms/QQ-plots and domain knowledge.

9) Optional: compare with SciPy

If you have SciPy available, you can verify that the NumPy implementation matches scipy.stats.skewtest.

x = rng.lognormal(mean=0.0, sigma=0.8, size=200)

ours = dagostino_skewtest(x)

try:
    from scipy.stats import skewtest

    scipy_res = skewtest(x)
    comparison = {
        "ours": ours,
        "scipy": {"z": float(scipy_res.statistic), "pvalue": float(scipy_res.pvalue)},
        "abs_diff_z": float(abs(ours["z"] - scipy_res.statistic)),
        "abs_diff_p": float(abs(ours["pvalue"] - scipy_res.pvalue)),
    }
except Exception as e:
    comparison = {"ours": ours, "scipy": None, "error": repr(e)}

comparison

{'ours': {'n': 200,
  'skewness': 2.5556762139133857,
  'z': 9.304601238442498,
  'pvalue': 0.0,
  'alternative': 'two-sided'},
 'scipy': {'z': 9.304601238442498, 'pvalue': 1.3449606122843268e-20},
 'abs_diff_z': 0.0,
 'abs_diff_p': 1.3449606122843268e-20}

Exercises

1. Simulate data from a symmetric but heavy-tailed distribution (e.g. Student t) and see how often the skewness test rejects.

2. Create a mixture of two normals with different means and observe how skewness changes.

3. For a fixed skewness, increase (n) and see how quickly the p-value shrinks.

References

• R. B. D’Agostino et al. (1990), A suggestion for using powerful and informative tests of normality

• SciPy documentation for scipy.stats.skewtest